I ran different models per trait with un-transformed data, log-transformed and sqrt-transformed data. Then I used a custom diagnostic function and the package DHARMa to compare between models. Then I selected the model that better fit the assumptions. The model I will show here are the ones selected. Finally, I used the package emmeans to estimate the marginal means per model.
This is the first model that looks at the differences between Tribulus mericarps, flowers and leaves from mainland and island populations.
## [1] "Kurtosis=0.819189880248705"
## [1] "Skew=-0.29428570229393"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(length)
## Chisq Df Pr(>Chisq)
## mainland_island 9.3113 1 0.002277 **
## year_collected 9.9210 1 0.001634 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=0.915288926843053"
## [1] "Skew=-0.184236416326029"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(width)
## Chisq Df Pr(>Chisq)
## mainland_island 10.36 1 0.001288 **
## year_collected 7.84 1 0.005110 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=0.989267427953835"
## [1] "Skew=-0.241253118215469"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: depth
## Chisq Df Pr(>Chisq)
## mainland_island 50.871 1 9.865e-13 ***
## year_collected 20.387 1 6.325e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=2.58419661904014"
## [1] "Skew=-0.421423446072704"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: spine_length
## Chisq Df Pr(>Chisq)
## mainland_island 0.3219 1 0.5705
## year_collected 24.7144 1 6.649e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=2.81329973280421"
## [1] "Skew=-0.769988792758626"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: tip_distance
## Chisq Df Pr(>Chisq)
## mainland_island 2.0743 1 0.1498
## year_collected 1.5347 1 0.2154
For spine number, I think it is a count trait, so I used this model:
glm(spine_num ~ mainland_island + year_collected, data = meri_spine.number, family = poisson)
But the diagnostic of this model looks different. How should I evaluate de assumptions of this model?
## [1] "Kurtosis=-0.0210478608607598"
## [1] "Skew=-0.395691016281737"
## Analysis of Deviance Table (Type II tests)
##
## Response: spine_num
## LR Chisq Df Pr(>Chisq)
## mainland_island 154.622 1 < 2.2e-16 ***
## year_collected 7.279 1 0.006975 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For lower spines, I think it is a count trait, so I used this model:
meri_lower.spines_m1 <- glm(lower_spines ~ mainland_island + year_collected, data = meri_lower.spines, family = “binomial”)
## [1] "Kurtosis=-1.61184279943993"
## [1] "Skew=-0.38863425255693"
## Analysis of Deviance Table (Type II tests)
##
## Response: lower_spines
## LR Chisq Df Pr(>Chisq)
## mainland_island 512.90 1 < 2.2e-16 ***
## year_collected 10.75 1 0.001042 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=1.42347292893909"
## [1] "Skew=0.00136618490346419"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: petal_length
## Chisq Df Pr(>Chisq)
## mainland_island 0.0000 1 0.99656
## year_collected 3.7983 1 0.05131 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=0.549228998714279"
## [1] "Skew=-0.231629387219308"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log(leaf_length)
## Chisq Df Pr(>Chisq)
## mainland_island 24.6162 1 6.996e-07 ***
## year_collected 5.3636 1 0.02056 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=0.704811597680861"
## [1] "Skew=-0.168061416832369"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log(leaflet_length)
## Chisq Df Pr(>Chisq)
## mainland_island 5.4526 1 0.019539 *
## year_collected 7.2378 1 0.007138 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For leaflet number, I think it is a count trait, so I used this model:
glm(number_of_leaflets ~ mainland_island + year_collected, family = poisson, data=leaf_length)
## [1] "Kurtosis=2.12466116567275"
## [1] "Skew=-0.651977029226767"
## DHARMa:plot used testOutliers with type = binomial for computational reasons (nObs > 500). Note that this method may not have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
## Analysis of Deviance Table (Type II tests)
##
## Response: number_of_leaflets
## LR Chisq Df Pr(>Chisq)
## mainland_island 33.467 1 7.25e-09 ***
## year_collected 1.189 1 0.2755
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## mainland_island response SE df asymp.LCL asymp.UCL
## island 6.06 0.0788 Inf 5.91 6.22
## mainland 5.71 0.0897 Inf 5.53 5.88
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
## mainland_island response SE df asymp.LCL asymp.UCL
## island 3.12 0.0369 Inf 3.05 3.20
## mainland 2.95 0.0426 Inf 2.86 3.03
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
## mainland_island emmean SE df asymp.LCL asymp.UCL
## island 4.75 0.0466 Inf 4.66 4.84
## mainland 4.25 0.0547 Inf 4.14 4.35
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## mainland_island emmean SE df asymp.LCL asymp.UCL
## island 4.30 0.106 Inf 4.09 4.51
## mainland 4.39 0.129 Inf 4.14 4.64
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## mainland_island emmean SE df lower.CL upper.CL
## island 7.93 0.202 277 7.53 8.33
## mainland 8.35 0.218 284 7.92 8.78
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## mainland_island rate SE df asymp.LCL asymp.UCL
## island 2.88 0.0401 Inf 2.80 2.96
## mainland 3.75 0.0593 Inf 3.63 3.87
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## mainland_island emmean SE df asymp.LCL asymp.UCL
## island 0.136 0.0326 Inf 0.0716 0.199
## mainland 2.086 0.0970 Inf 1.8962 2.277
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
In general mericarps from mainland populations are smaller than island mericarps. However, spine size seems to differ with mainland populations having longer tip distances and higher spine numbers than island populations.
## mainland_island emmean SE df lower.CL upper.CL
## island 16.8 0.285 396 16.2 17.3
## mainland 16.8 0.308 377 16.2 17.4
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## mainland_island response SE df lower.CL upper.CL
## island 29.3 0.765 384 27.8 30.8
## mainland 24.0 0.760 417 22.5 25.5
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## mainland_island response SE df lower.CL upper.CL
## island 8.29 0.179 391 7.94 8.64
## mainland 7.67 0.200 417 7.29 8.07
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## mainland_island emmean SE df asymp.LCL asymp.UCL
## island 2.66 0.0104 Inf 2.64 2.68
## mainland 2.55 0.0154 Inf 2.52 2.58
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
Mainland populations seems to have smaller leaves and less leaflets than island populations.
This is the second model, that shows the differences of Tribulus traits from Galapagos populations and other Islands.
## [1] "Kurtosis=1.82233135992289"
## [1] "Skew=-0.357589730850339"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(petal_length)
## Chisq Df Pr(>Chisq)
## galapagos_other 87.9227 1 < 2.2e-16 ***
## year_collected 8.8887 1 0.002869 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=0.174591449444752"
## [1] "Skew=-0.094965010791886"
## qu = 0.75, log(sigma) = -2.141644 : outer Newton did not converge fully.
## qu = 0.75, log(sigma) = -2.176578 : outer Newton did not converge fully.
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log(leaf_length)
## Chisq Df Pr(>Chisq)
## galapagos_other 18.0889 1 2.108e-05 ***
## year_collected 0.5641 1 0.4526
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=0.295097534906618"
## [1] "Skew=-0.0991019439431268"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log(leaflet_length)
## Chisq Df Pr(>Chisq)
## galapagos_other 17.3782 1 3.063e-05 ***
## year_collected 1.0138 1 0.314
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For leaflet number, I think it is a count trait, so I used this model:
glm(number_of_leaflets ~ mainland_island + year_collected, family = poisson, data=leaf_length)
## [1] "Kurtosis=1.11565272013564"
## [1] "Skew=-0.378042815760105"
## DHARMa:plot used testOutliers with type = binomial for computational reasons (nObs > 500). Note that this method may not have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
## Analysis of Deviance Table (Type II tests)
##
## Response: number_of_leaflets
## LR Chisq Df Pr(>Chisq)
## galapagos_other 3.3502 1 0.06720 .
## year_collected 3.3740 1 0.06623 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## galapagos_other response SE df lower.CL upper.CL
## Galapagos 9.03 0.705 219 7.7 10.5
## other 17.15 0.269 211 16.6 17.7
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
Petal length from Galapagos seems to be smaller than other island systems.
## galapagos_other response SE df lower.CL upper.CL
## Galapagos 24.1 1.155 158 21.9 26.5
## other 30.9 0.945 265 29.1 32.8
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## galapagos_other response SE df lower.CL upper.CL
## Galapagos 7.02 0.293 169 6.47 7.63
## other 8.66 0.226 260 8.23 9.12
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## galapagos_other emmean SE df asymp.LCL asymp.UCL
## Galapagos 2.64 0.0172 Inf 2.60 2.67
## other 2.68 0.0160 Inf 2.65 2.71
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
Leaf, leaflet are smaller in the Galapagos compared to other islands. Leaflet number are also fewer in Galapagos.
For this third model I filter the samples from Galapagos only. Then I defined the Islands that have different finch communities based on the presence or absence of large ground finches: G. magnirostris, G. cornirostris.
Floreana, San Cristobal, Santa Fe, Champion, Baltra, Enderby, Gardner and Daphne Major, previous 1983 are considered without magnis or cornirostris.
The rest of the Islands are considered with large ground finches.
## [1] "Kurtosis=0.640272917671844"
## [1] "Skew=0.0069518899108278"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: length
## Chisq Df Pr(>Chisq)
## finch_beak 5.0459 1 0.02468 *
## year_collected 0.0972 1 0.75522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=0.866086533593143"
## [1] "Skew=-0.13190943721093"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(width)
## Chisq Df Pr(>Chisq)
## finch_beak 0.5813 1 0.4458
## year_collected 1.7970 1 0.1801
## [1] "Kurtosis=0.741297164859156"
## [1] "Skew=-0.20439673252663"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: depth
## Chisq Df Pr(>Chisq)
## finch_beak 0.2483 1 0.61824
## year_collected 3.9913 1 0.04574 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=2.22346563111609"
## [1] "Skew=-0.405909926348517"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: spine_length
## Chisq Df Pr(>Chisq)
## finch_beak 0.0991 1 0.7529296
## year_collected 14.3737 1 0.0001499 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=2.76047469083229"
## [1] "Skew=-0.776310384626465"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: tip_distance
## Chisq Df Pr(>Chisq)
## mainland_island 2.0743 1 0.1498
## year_collected 1.5347 1 0.2154
For spine number, I think it is a count trait, so I used this model:
glm(spine_num ~ mainland_island + year_collected, data = meri_spine.number, family = poisson)
But the diagnostic of this model looks different. How should I evaluate de assumptions of this model?
## [1] "Kurtosis=-0.622289864714663"
## [1] "Skew=0.211561144749239"
## DHARMa:plot used testOutliers with type = binomial for computational reasons (nObs > 500). Note that this method may not have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
## Analysis of Deviance Table (Type II tests)
##
## Response: spine_num
## LR Chisq Df Pr(>Chisq)
## finch_beak 0.0408 1 0.83990
## year_collected 6.5550 1 0.01046 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For lower spines, I think it is a count trait, so I used this model:
meri_lower.spines_m1 <- glm(lower_spines ~ mainland_island + year_collected, data = meri_lower.spines, family = “binomial”)
## [1] "Kurtosis=-1.86746698858495"
## [1] "Skew=-0.0114759981319253"
## Analysis of Deviance Table (Type II tests)
##
## Response: lower_spines
## LR Chisq Df Pr(>Chisq)
## finch_beak 193.586 1 < 2.2e-16 ***
## year_collected 11.904 1 0.0005601 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=-0.98535796283949"
## [1] "Skew=0.165132634915697"
## qu = 0.25, log(sigma) = -3.886972 : outer Newton did not converge fully.
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: petal_length
## Chisq Df Pr(>Chisq)
## finch_beak 0.5165 1 0.4724
## year_collected 6.2386 1 0.0125 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=-0.0288924193805786"
## [1] "Skew=0.340587284158791"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(leaf_length)
## Chisq Df Pr(>Chisq)
## finch_beak 1.2475 1 0.2640318
## year_collected 13.6747 1 0.0002174 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Kurtosis=-0.0292060182499396"
## [1] "Skew=-0.135277570705385"
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log(leaflet_length)
## Chisq Df Pr(>Chisq)
## finch_beak 1.6469 1 0.19938
## year_collected 4.8053 1 0.02837 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For leaflet number, I think it is a count trait, so I used this model:
glm(number_of_leaflets ~ mainland_island + year_collected, family = poisson, data=leaf_length)
## [1] "Kurtosis=1.0223498497344"
## [1] "Skew=0.114640438164658"
## Analysis of Deviance Table (Type II tests)
##
## Response: number_of_leaflets
## LR Chisq Df Pr(>Chisq)
## finch_beak 0.01639 1 0.8981
## year_collected 2.17625 1 0.1402
## finch_beak emmean SE df asymp.LCL asymp.UCL
## 0 5.98 0.146 Inf 5.69 6.27
## 1 6.41 0.131 Inf 6.15 6.67
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## finch_beak response SE df asymp.LCL asymp.UCL
## 0 3.10 0.0529 Inf 2.99 3.20
## 1 3.15 0.0487 Inf 3.06 3.25
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
## finch_beak emmean SE df asymp.LCL asymp.UCL
## 0 4.82 0.0744 Inf 4.67 4.96
## 1 4.77 0.0663 Inf 4.64 4.90
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## finch_beak emmean SE df asymp.LCL asymp.UCL
## 0 4.19 0.216 Inf 3.76 4.61
## 1 4.28 0.194 Inf 3.90 4.66
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## finch_beak emmean SE df lower.CL upper.CL
## 0 7.34 0.439 122 6.47 8.2
## 1 7.93 0.391 126 7.15 8.7
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## finch_beak emmean SE df asymp.LCL asymp.UCL
## 0 0.990 0.0215 Inf 0.948 1.03
## 1 0.996 0.0224 Inf 0.952 1.04
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## finch_beak emmean SE df asymp.LCL asymp.UCL
## 0 -0.479 0.0500 Inf -0.577 -0.381
## 1 0.474 0.0482 Inf 0.380 0.569
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
Within the Galapagos Islands Tribulus mericarps have differences associated to finch communities. Mericarps that are larger in length (more seeds) and have lower spines are associated with large finch communities.
Width, Depth, Spine Length, Tip Distance and Spine Number do not show a large difference between communities.
## finch_beak emmean SE df lower.CL upper.CL
## 0 10.29 1.27 20.4 7.65 12.9
## 1 9.22 0.86 20.3 7.43 11.0
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## finch_beak response SE df lower.CL upper.CL
## 0 21.8 2.00 47.4 18.0 26.0
## 1 24.4 1.29 40.5 21.8 27.1
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
## finch_beak response SE df lower.CL upper.CL
## 0 6.21 0.489 48.0 5.3 7.27
## 1 6.96 0.339 41.4 6.3 7.67
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## finch_beak emmean SE df asymp.LCL asymp.UCL
## 0 2.65 0.0302 Inf 2.59 2.71
## 1 2.65 0.0171 Inf 2.62 2.69
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95